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ABSTRACT 

Most ultraluminous X-ray sources (ULXs) are believed to be X-ray binary 
systems, but previous observational and theoretical studies tend to prefer a black 
hole rather than a neutron star accretor. The recent discovery of 1.37 s pulsations 
from the ULX M82 X-2 has established its nature as a magnetized neutron star. 

In this work we model the formation history of neutron star ULXs in an M82- or 
Milky Way-like galaxy, by use of both binary population synthesis and detailed 
binary evolution calculations. We hnd that the birthrate is around 10“^ yr“^ for 
the incipient X-ray binaries in both cases. We demonstrate the distribution of 
the ULX population in the donor mass - orbital period plane. Our results suggest 
that, compared with black hole X-ray binaries, neutron star X-ray binaries may 
significantly contribute to the ULX population, and high-mass and intermediate- 
mass X-ray binaries dominate the neutron star ULX population in M82- and 
Milky Way-like galaxies, respectively. 

Subject headings: accretion - stars: neutron - stars: evolution - X-rays: binaries 


1. Introduction 


Ultraluminous X-ray sources (ULXs) are off-nuclear, point-like sources with X-ray lumi 
nosities exceeding 10^® ergs“^, first discovered in nearby galaxies with Einstein (Fabbiano 


1989). More recent observations with improved X-ray telescopes, such as Chandra and 


XMM-Newton, have greatly increased the number of this kind of sources (Fabbiano & White 


2006; Roberts 2007; Feng & Soria 2011, for reviews). They are most likely X-ray binaries 
(XRBs), in which a compact object accretes from a donor star through Roche-lobe over¬ 
flow (RLOF), but the nature of these objects has not been completely uncovered. If the 
radiation is isotropic and below the Eddington limit, the extremely high luminosities imply 
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the presence of an accreting intermediate-mass (10^ — 10^ M©) black hole (BH) (Colbert 


& Mushotzky 1999). Alternatively they are believed to be stellar-mass BHs with super- 


Eddington accretion. For example, King et al. (2001) proposed that geometrical beaming in 


the case of rapid accretion could lead to a very high apparent luminosity for a stellar-mass 


BH. Begelman (2002) showed that the isotropic luminosity of an accreting BH can exceed 


the Eddington limit by a factor of about 10 due to the photon-bubble instability in the 
accretion disk. Recently, the masses of two ULXs were dynamically measured to be in the 


stellar-mass BH range (Liu et al. 2013; Motch et al. 2014). 


Since the X-ray luminosities of ULXs are significantly higher than the Eddington limit 
Le (around 2 x 10^®erg s“^) for a 1.4M0 neutron star (NS(Q it is challenging to explain these 
bright X-ray sources with an accreting NS. However, the discovery of 1.37 s pulsations from 


the ULX M82 X-2 has provided unambiguous evidence for its NS nature (Bachetti et al. 


2014), implying that accreting NSs also contribute to the ULX population. 


M82 X-2 is in a 2.53 day orbit with a companion star more massive than 5.2 Mq 
( B achetti et al.||2014 ). The extremely high (isotropic) luminosity (around 10^° ergs“^) (Feng 

(Bachetti et al. 2014) clearly 


et al. 


2010) and rapid spin-up (at a rate of —2 x 10 


-10 


ss 


indicate that the binary is undergoing rapid mass transfer. Since the companion star is 
significantly more massive than the NS, the mass transfer is subject to delayed dynamical 
instability, and must be currently in the early phase of RLOF. After that the NS will be 
engulfed by the transferred matter, resulting in the formation of a common envelope (CE) 


(Bhattacharya & van den Heuvel 1991). The objective of this paper is to investigate how 


many such accreting NSs can be responsible for ULXs in a galaxy like M82 or the Milky Way 
(MW). We first explore the properties of the incipient NS XRBs, using a binary population 
synthesis (BPS) method (in Sect. 2), then calculate the detailed evolutions of these XRBs to 
obtain the numbers and the luminosity functions of the ULXs (in Sect. 3). We summarize 
our results in Sect. 4. 


2. Generation of the incipient NS XRBs 


To model the formation history of NS XRB-ULXs, we adopt the BPS code initially 
developed by Hurley et al. (2002) to calculate the evolution of a large population of the 


primordial binaries. We have updated and modified the code in several aspects (see Shao & 


^^Jfthe^Jfi_£ossesses_astrongjnagneti£fi£ldtodmnneHhe^ccretingjnaterialmntoJts_^urfecej_rtie_critic^ 

luminosity could be higher, about ~ 4(k^)I,g, where Iq and do are the length and the thickness of 

the accreting funnel, respectively ([Basko &: Sunyaev||1976p. 
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Li 2014, for details), especially the conditions for dynamically stable mass transfer and the 


treatments of CE evolution, which are briefly described as follows. 


During the evolution of a primordial binary, the primary first evolves off the main 
sequence and expands in size. This can lead to RLOF onto the secondary, causing it to be 
spun up and rejuvenated. If the mass transfer proceeds so slowly that the secondary can 
remain in thermal equilibrium, the mass transfer is thought to be stable. Otherwise the 
secondary will get out of thermal equilibrium and expand. This expansion may Anally cause 


the secondary to All its own RL, leading to the formation of a contact binary (Nelson &: 


Eggleton 2001). A critical mass ratio gcr is usually used to determine whether or not the 


mass transfer is dynamically stable in a binary. Instead of using the empirical results for 
gcr of Hurley et al. (2002), Shao & Li (2014) numerically calculate it considering both the 


response of the secondary to mass accretion and the effect of possible mass loss. Here we 


use the results of Model H given by Shao & Li (2014), in which it is assumed that half of 


the transferred mass is accreted by the secondary, and the other half is lost from the system. 
This model can well fit the observational distribution of Be/X-ray binaries. 

When the mass transfer is dynamically unstable, we employ the standard energy con¬ 


servation equation (Webbink 1984) to deal with the subsequent CE evolution, that is, in 


the spiral-in phase, the orbital energy of the secondary is used to expel the envelope of the 
primary. When calculating the binding energy of the stellar envelope, we include the con¬ 


tribution from the internal energy and adopt the fitting formulae of Xu & Li (2010) for the 


binding energy parameter A. We assume the CE efficiency ctce = 1-0 in our calculations. 
The evolution of a binary is determined by the primary mass Mi, secondary mass M 2 , 


and orbital angular momentum (Hurley et al. 2002). The orbit is invariably circularized 


before interaction by standard tidal interactions, so all binary orbits are taken to be circular. 


We adopt the initial mass function given by Kroupa et al. (1993) for the primary stars, and 


a flat distribution between 0 and 1 for the initial mass ratio of the secondary to the primary. 
For the initial orbital separation a, we assume that In a is evenly distributed between a = 3 Rq 
and 10 ^Rq. We adopt solar metallicity Z = 0.02 in our simulations. 


We consider two cases of star formation activities for star forming galaxies like M82 and 
late-type galaxies like the MW. We adopt a constant overall star formation rate of IOMq yr“^ 
over the last 100 Myr (Forster Schreiber et al. 2003; Yao 2009[ ) in case (1) , and a constant 
star formation rate 3MQyr“^ over the 13 Gyr period in case (2). 


We have evolved a population of 7 x 10® primordial binaries, and generated a subset of 
about 1.0 X 10^ and 1.5 x 10^ incipient XRBs containing a NS and an unevolved secondary 
star of mass lower than 20Mq in cases (1) and (2), respectively. Note that for the mechanisms 
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of NS formation we consider both core-collapse supernovae of massive stars and electron- 
capture supernovae of intermediate-mass stars, following the criterion suggested by |Fryer| 
et al.| ( |2012| ). In Fig. 1 we plot the distribution of these XRBs in the donor (or secondary) 
mass (Md) - orbital period (Porb) plane (the left panel) and the birthrate distribution as a 
function of the donor mass and of the orbital period (right panel). Note that in case (1) 
most of the secondaries are more massive than SMq (upper panel), because of the much 
shorter evolution time than in case (2) (lower panel). The total birthrates are estimated to 
be about 1.6 x 10“^ yr“^ and 6.9 x 10“^ yr~\ respectively. 


3. Evolution to ULXs 


Based on the result in Fig. 1 we calculate the evolution of the generated NS XRBs with 


the TWIN version of the stellar evolution code developed by Eggleton (1971, 1972). Here 


the initial mass of the NS is assumed to be 1.4M0. We have evolved thousands of binary 
systems with the donor mass varying from O.SMq to 20Mq by steps of 0.25Mq, and the 
orbital period Porb (in units of days) increasing logarithmically from —0.5 to 3 by steps of 
0.1, as obtained from Fig. 1. We use these binaries to represent all the XRBs, and their 
numbers are calculated by cumulating the XRBs in a specific matrix of AM^ x A (log Porb) 
in the — Porb plane by weighing their formation rates and life spans. Note that the XRBs 
are initially eccentric due to the supernova kicks, which are assumed to have a Maxwellian 
distribution with a dispersion a — 265 kms“^ ( Hobbs et al.|2005 ) for core-collapse supernovae 
and 50 kms“^ (Dessart et ah 2006) for electron-capture supernovae. Here we assume that the 


orbital angular momentum of an incipient XRB is conserved and then the system is quickly 
circularized with a new separation, which is smaller by a factor of (1 — e^) (Bhattacharya & 


van den Heuvel 1991) 


In Eggleton’s code mass transfer via RLOE can be modeled as a function of the potential 
difference A(/) between the stellar surface and the RL surface. When one of the stars overfills 
its RL, the mass flux at each mesh point outside the Roche surface is given by 


dM 

dm 


= - 10 ' 




( 1 ) 


where m and r are the mass coordinate and the radius, respectively. The mass transfer rate 
is then calculated by integrating Eq. (1) over all mesh points outside the Roche surface 
potential. Actually before the photospheric radius of a massive donor star reaches its RL, 
the atmospheric matter begins to spill over towards the NS along the inner Langrangian 
point. This phase of beginning atmospheric RLOE precedes the main phase of RLOE until 


the mass transfer rate rises to the Eddington value (Savonije 1979). Subsequently the mass 
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transfer rate increases rapidly to become super-Eddington. In Fig. 2 we show the exampled 
evolutionary tracks of three binary systems. The initial parameters are = 6 Mq and 
-Porb = 1 d, Md = 6Mq and Porb = 10 d, and = IOMq and Porb = 1 d in the top, middle, 
and bottom panels, respectively. In the top panel, the donor evolves to overflow from the 
RL at the age of 39.4 Myr. The mass transfer rate increases from 2.1 x lO“^°M0yr“^ to 
2.4 X lO“^M0yr“^ within a time of 4.46 x 10^ yr. The orbital period decays to 0.58 d, and 
the donor mass decreases to 5.59Mq, which means that O.4IM0 of the donor’s envelope is 
stripped before the CE occurrence. In the middle panel, a longer initial orbital period of 10 
d is set for the binary. The onset of RLOF occurs at the age of 66.08 Myr, and the donor star 
is more evolved. The mass transfer rate rises from 1.3 x 1O“^°M0 yr“^ to 1.1 x 1O“^M0 yr“^ 
within a shorter time of 1.2 x 10^ yr. The orbital period drops to 3.68 d, and 0.78Mq material 
is transferred from the donor within this time. In the bottom panel, the initial donor star is 
more massive (with a mass of IOM0), and the mass transfer lasts 2.32 x 10^ yr. The orbital 
period decreases to 0.75 d when 0.17 Mq of the donor’s envelope is transferred. 

Given the mass transfer rate, we calculate the X-ray luminosity in two ways. In the first 
one, we calculate it with the traditional formula 

Lx = O.lMc^ (2) 


without considering the Eddington limit. In the second, we use the same formula for sub- 
Eddington accretion rates. When M is higher than the Eddington accretion rate Mg, we 


adopt the model of King (2008) to convert the mass transfer rates into the X-ray luminosities. 


In this model, the accretion disk becomes geometrically thick, which influences the X-ray 
luminosity in two ways. First, radiation becomes less efficient and the bolometric luminosity 
no longer follows M linearly. Second, the outgoing radiation may be collimated due to a 
biconical geometry at the inner part of the accretion disk. The accretion luminosity is then 
contributed by two parts. The region outside the so-called spherization radius Pgph where 
the mass inflow first becomes locally Eddington ( Shakura fc Syunyaev||1973 ; Begelman et al. 
2006) releases the accretion luminosity close to Lg. For the region within Pgph, the accretion 


luminosity is about ln(Psph/3Ps) ~ ln(M/ME), where Rs is the Schwarzschild radius (Frank 


et al. 2002). The total luminosity is then (King 2008), 


L« 


Lf 


1 -|- In 


M 


( 3 ) 


Because of the geometric collimation, one can see the source in directions within one of the 
cones, with an apparent (isotropic) X-ray luminosity 


T 

Lx^ -r- 


1 -|- In 


M 


( 4 ) 
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where b is the beaming factor, possibly depending on M. Here we adopt 6 ~ 0.1 as suggested 
by|Kmi(p008l). 


In Fig. 3 we present the number distribution of the predicted ULX population in case 
(1) (i.e., in a galaxy like M82) as a function of and Porb, with Lx greater than 10^® ergs“^. 
In the upper and lower panels the X-ray luminosities are calculated with Eqs. (2) and (4), 
respectively. As noted before, for each matrix element in the left panel, the number is 
calculated by multiplying the birthrate of the incipient XRBs with the evolutionary time 
span within the matrix element. It is seen that most ULXs tend to be high-mass XRBs 
in relatively short orbits (with orbital periods shorter than a few days). Figure 4 shows 
the same distributions in case (2) (i.e., in a MW-like galaxy). A comparison of Figs. 3 and 
4 shows that in the latter case XRBs with donor of mass lower than 3Mq dominate the 
population because of their much longer lifetime. 


In Fig. 5 we plot the X-ray luminosity function of the ULX population in cases (1) (left) 
and (2) (right). The red and black curves correspond to higher than 2Mq and 5Mq, 
respectively. For the dashed and solid curves the luminosities are calculated with Eqs. (2) 
and (4), respectively. In case (2) the effect of anisotropic radiation on the observable number 
is taken into account. The predicted ULX numbers lie between a few tenths and a few in 
each case. Madhusudhan et ah (2008) investigated the evolution of ULXs consisting of a 
stellar-mass BH accretor, and found that their numbers range from 0.1 to 0.2 for Lx higher 
than 2 x 10^® ergs“^ in a MW-like galaxy. Hence accreting NSs may play an important role 
in the formation of ULXs, compared with BHs. 


4. Summary 

We have calculated the formation history of NS ULXs. Our results show that, (1) 
compared with BH XRBs, NS XRBs may significantly contribute to the ULX population; 
(2) high-mass and intermediate-mass XRBs dominate the ULX population in M82- and 
MW-like galaxies, respectively. 


This work was supported by the Natural Science Eoundation of China under grant 
numbers 11133001, 11203009, and 11333004, the Strategic Priority Research Program of 
CAS (under grant number XDB09000000), and the graduate innovative project of Jiangsu 
Province (CXZZ13-0043). 
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Fig. 1.— The upper and lower panels show the predicted distributions of the incipient NS 
XRBs and their birthrates in cases (1) and (2), respectively. The left panel shows the XRB 
distribution in the donor mass vs. orbital period plane. Their birthrate distributions are 
presented in the right panel. 
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Fig. 2.— Example evolution of the orbital period and the mass transfer rate for three NS 
XRBs, as a function of the donor mass and the age, respectively. The initial parameters are 
Md = 6M0 and Porb = 1 d, Md = QMq and Porb = 10 d, and Md = IOMq and Porb = 1 d in 
the top, middle, and bottom panels, respectively. 
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Fig. 3.— The predicted distributions of NS ULXs in case (1). In the upper and lower 
panels, the X-ray luminosities are calculated with Eqs. (2) and (4), respectively. The left 
panel shows the ULX distribution in the — Porb plane. The color in each matrix element 
represents the number of ULXs with Lx higher than 10^® ergs“^. The right panel shows the 
number distribution of the ULXs as a function of the donor mass and the orbital period 
Porb • 
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Fig. 4.— Same as Fig. 3, but for case (2). 
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Fig. 5.— The left and right panels show the X-ray luminosity functions of the NS ULX 
population in case (1) and (2), respectively. The red and black curves correspond to the 
donor mass higher than 2Mq and 5Mq, respectively, the X-ray luminosities for the dashed 
and solid curves are calculated with Eqs. (2) and (4), respectively. 







